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CNl ' Abstract 
>■ 

' We present several topics involving the computation of dynamical systems. The em- 

phasis is on work in progress and the presentation is informal - there are many technical 
Cp ■ details which are not fully discussed. The topics are chosen to demonstrate the various 

interactions between numerical computation and mathematical theory in the area of 
dynamical systems. We present an algorithm for the computation of stable manifolds 
Q>^ [ of equilibrium points, describe the computation of Hopf bifurcations for equilibria in 

■ parametrized families of vector fields, survey the results of studies of codimension two 
global bifurcations, discuss a numerical analysis of the Hodgkin and Huxley equations, 

■ and describe some of the effects of symmetry on local bifurcation. 
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^ ■ 1 Introduction 



This article is a written record of the lectures delivered at the 1992 NATO Summer School 
in Montreal. The spirit of the lectures was informal, and we have tried to preserve the 
informality and didactic quality of the lectures. The lectures dealt with several related 
I topics, all involving computations of dynamical systems. The emphasis here stresses topics 

that have not yet been fully developed in other recent papers. Accordingly, they reflect 
to a large extent "work in progress" rather than presenting a summary and review that 
makes a finished mathematical tale. There is also a lack of rigor that reflects a pragmatic 
approach that seeks reliable answers to questions even when these answers cannot be totally 
integrated into formal mathematical theory. 

Our goal is to answer questions about speciflc dynamical systems. If / : i?" x — > i?" 
defines a k parameter family of vector fields, we would like to know everything about the 
phase portraits of / and the bifurcations that occur as parameters are varied. Since most 
dynamical systems cannot be integrated in terms of analytic expressions that are valid for 
all time, this is a task that can only be approached sensibly through numerical computa- 
tion. This is problematic for mathematicians because naive error estimates tend to grow 
exponentially with iterative calculations such as the numerical integration of a vector field. 
For the most part, we ignore this difficulty and proceed with confidence that computers 
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provide good approximations to trajectories of vector fields. We seek to find additional al- 
gorithms that allow us to calculate aspects of dynamical systems that are hard to determine 
otherwise. Though these lectures do not emphasize applications, our interests extend be- 
yond individual algorithms to the implementation of efficient and powerful problem solving 
environments for the analysis of dynamical systems and their bifurcations |^. 

2 Computing Stable Manifolds of Equilibria 

The stable and unstable manifolds of equilibrium points are important geometric objects in 
the phase portraits of vector fields. The Stable Manifold Theorem asserts their existence, 
and some proofs of the Stable Manifold Theorem are constructive; e.g. [^3|| . Nonethe- 
less, the numerical computation of stable and unstable manifolds of dimension larger than 
one presents difficulties. This section describes some of these difficulties and sketches an 
approach that leads to a reasonable algorithm for displaying the two-dimensional stable 
manifold of the Lorenz system l50| ]. 

We begin by recalling the Stable Manifold Theorem for equilibria of finite-dimensional 
vector fields: 

Theorem 2.1 Let X he a vector field on an n- dimensional manifold M with an equilib- 
rium point p and flow <I>. Let L : Tp{M) Tp{M) he the linearization of X at p. Assume 
that Tp{M) splits as an invariant direct sum + with the spectrum of L\E^ in the 
open left half plane of C and the spectrum of L\E^ in the open right half plane of C. Then 
there are C' suhmanifolds W and passing through p with tangent spaces E" and E^ 
respectively with the following property: 

W" = {x £ M\<^>{t, x) ^p as t^ oo} 

W"" = {x £ M|$(t, x)^past^ -oo} . 

If ^ is complete, then and are 1 — 1 immersions of Euclidean spaces into M. 

There is a naive approach to computing the stable manifold of an equilibrium that can 
be derived from the simplest proofs of the Stable Manifold Theorem. If p is an equilibrium 
point with the linearization at p having stable manifold E'^ and unstable manifold the 
stable manifold of p can be represented in local coordinates as the graph of a function 
7 : E^** — > S". Near p, each affine subspace parallel to intersects in a single point. 
The function 7 has vanishing derivative at p, so points of E^ near p form a good 
approximation to there. Using the invariance of the stable manifold, the backwards 
trajectories starting on a sphere surrounding p in E^ sweep out an approximation to W^. 
One can choose a set of initial conditions on E'^ and compute their trajectories to visualize 
ly*, but the process often works poorly. 

Let us examine the Lorenz system as an example. The equations of the Lorenz system 

are 

X = a{y — x) 
y = px — y — xz 
z = —/3z + xy 
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with "standard" parameters a = 10, /? = 8/3, and p = 28. For these parameters, the origin 
is a saddle point with a two-dimensional stable manifold and a one-dimensional unstable 
manifold. The linearization at the origin is given by the matrix 

/ -10 10 \ 

28-1 
^ -8/3 / 

The negative eigenvalues of this matrix are approximately —2.67 and —22.8 with a ratio 
that is approximately 8.56. There are two problems associated with the large ratio between 
these eigenvalues. First, "most" trajectories in the stable manifold approach the origin from 
directions close to the z axis which is the eigendirection of the eigenvalue —8/3. This can 
be seen readily by solving the linear system 

X = —ax 
y = -by . 

The trajectories of this system lie along the curves y = cx""!^ . Points uniformly distributed 
on the unit circle approach the origin with much higher density close to the y axis if a <C 6. 
The second problem is that the trajectories flow much faster in the eigendirection of the 
eigenvalue of larger magnitude than in the eigendirection of the weaker eigenvalue. A small 
circle of initial conditions in the stable manifold will flow backwards to an ellipse whose axes 
have a ratio approximately e*^'^^ ~ 5220 after one unit of time. Thus a fixed time integration 
of the vector field on a set of initial conditions clustered near the weak eigendirection still 
will yield a set of points that stretches along the strong eigendirection. 

This discussion indicates that strategies for drawing stable manifolds based solely upon 
the numerical integration of trajectories starting from initial values near the equilibrium 
may be ineffective. To obtain approximations to large regions of a stable manifold, our 
approach is to mimic the construction of "geodesic coordinates" in differential geometry. 
We describe the general construction. 

Let p be a hyperbolic equilibrium point of a vector field X on a Riemannian manifold 
M and let W = W^{p) be the stable manifold of p. We make use of the Riemannian 
metric induced on W from its embedding in the ambient Riemannian manifold M. There 
is a neighborhood U oi p in W such that U — {p} is foliated by spheres Sr which lie at a 
distance r from p in the Riemannian metric. We would like to write a differential equation 
for the evolution of these spheres with increasing r, and this is easy to do when the vector 
field X is not tangent to Sr- The fundamental observation that allows us to proceed is that 
the tangent space Tx{W) to W at x G Sr is spanned by X and Tx{Sr) as long as the vector 
field is not tangent to Sr- The unique geodesic in W from p to any point of the sphere Sr is 
characterized by the requirement that it is orthogonal to Sr and that its tangent has unit 
length and lies in the subspace spanned by the tangents to Sr and the vector field X- This 
provides the motivation for an algorithm to evolve the sphere Sr by extending the geodesies 
through p and points on the sphere, as pictured in Figure ||. 

Theorem 2.2 Let X be a C^' vector field on an n- dimensional Riemannian manifold M 
with a hyperbolic equilibrium point p and flow Let W be the s-dimensional stable manifold 
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of p and let Sr be the geodesic sphere at distance r from p in the Riemannian metric induced 
on W. If X is not tangent to Sr, then the vector field obtained by orthogonally projecting 
X onto the normal bundle of Sr is tangent to the geodesic rays of W emanating from p. 

There may well be tangencies between the vector field X and the geodesic spheres Sr in 
a stable manifold. When this happens, it is no longer a simple matter to "evolve" Sr as a 
submanifold without additional geometric information. We approach the problem with the 
question of determining when does the set of trajectories passing through a submanifold 
S of initial conditions form a smooth submanifold W. We point out two difficulties. The 
following example illustrates one problem that occurs if there are isolated points of tangency 
between X and S. 

Let X be the vector field X = dx — xdy + x'^dz in and let S be the x axis. The vector 
field X is tangent to S at the origin. The flow of X is defined by the map ^{x,y,z,t) = 
{x + t,y — xt — z + x'^t + xt'^ + The surface of trajectories with initial conditions 
on the curve 5 is given by the map F{x,t) = {x + t, —xt — t'^/2,x'^t + xt'^ + t^/3). The image 
of F is a singular surface whose intersection with the plane x = is a cusp. This example 
indicates that given a smooth submanifold S of dimension s — 1, there are simple vector 
fields with the property that the trajectories passing through the submanifold do not form 
a smooth surface of dimension s. Locally, a submanifold S may be evolved by a vector field 
X in such a fashion, as long as the tangent spaces to S do not contain X. However, at 
points of tangency additional information is required to ensure that the trajectories through 
S will form a smooth submanifold. 

There is a further difficulty that we encounter in evolving submanifolds. Suppose that W 
is a submanifold with boundary and that X is a vector field with a segment of a trajectory 
7 in the boundary of W. Assume that W is contained in the interior of an invariant 
manifold W of the same dimension. It is not true that W and the vector field define W. If 
the trajectory 7 never enters the interior of W, then W — W may consist of points whose 
trajectories do not intersect W at all. For example, the most dramatic case occurs when the 
dimension of W is two and the boundary of is a periodic orbit. In this case, we cannot 
"grow" W in the fashion that we wish to use for stable and unstable manifolds. On the other 
hand, the scenario we have just described cannot occur for a stable or unstable manifold, 
since these invariant manifolds are formed by a set of trajectories that are asymptotic to 
an equilibrium point. 

This story ends in an unfinished state. We would like to use algorithms based upon 
computing geodesies for computing stable and unstable manifolds because they display the 
geometry of the manifolds in a way that naive computation of trajectories does not. On the 
other hand, there are mathematical difficulties with the formulation of algorithms that are 
based solely upon the computation of a vector field that will be tangent to geodesies. This 
is an area of research in which more work is warranted. The goal is to find a procedure 
that "works" for the largest possible class of examples, where the criterion of success is 
the ability to carry out the computation to a specified precision in a reasonable amount of 
time. Figure ^ shows a picture of the computed stable manifold for the origin of the Lorenz 
system. 
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3 Computing Hopf Bifurcations 

This section describes joint work in progress with Mark Myers and Bernd Sturmfels. The 
question we discuss involves the computation of points where Hopf bifurcations occur in 
parametrized famihes of vector fields. Let x = f\{x) be a system of differential equations 
on i?" depending upon a /c-dimensional parameter A. One of the first steps in a bifurcation 
analysis of this family is the determination of parameter values at which there are non- 
hyperbolic equilibrium points (x, A) G i?" x . These points are characterized by the 
conditions that f\{x) = and Df\[x) has a zero or purely imaginary eigenvalue. In 
both cases, the set of parameter values at which bifurcation occurs is expected to have 
codimension one in the parameter space. Generically, the case of a zero eigenvalue is a 
saddle-node bifurcation while the case of pure imaginary eigenvalues is a Hopf bifurcation 
p6| . Note that solving the equations describing either type of bifurcation does not involve 
numerical integration of the differential equations. Instead, one has a root-finding problem 
to solve. Since the Jacobian Df\{x) depends upon /, it is not immediately apparent that 
one can express this problem as a non-degenerate system of n-|-l equations 'm.n+k variables. 

For saddle-node bifurcations, zero eigenvalues of the Jacobian can be detected by com- 
puting the determinant of the Jacobian. It is not difficult to prove that the map F : 
j^n ^ ^fc _^ given by F{x,X) = {fx{x),det{Dfx{x)) is non-degenerate for generic / at 

most points [^]. The fundamental example is given by the normal form for a saddle-node 
bifurcation: f\{x) = A -|- x^. In this example, F{x, A) = (A -|- x^, 2x) is non-singular. The 
general saddle-node bifurcation can be transformed to this normal form by using the tech- 
niques of singularity theory [^2|. In implementing algorithms for computing saddle-node 
bifurcations, we face questions about the efficiency and accuracy with which the determinant 
of the Jacobian of a map can be computed. 

We seek algorithms that will locate Hopf bifurcations in a manner analogous to the 
computation of saddle-node bifurcations. Thus, we want a map H : x R^ — > 
which vanishes when Df has an equilibrium at which there are a pair of pure imaginary 
eigenvalues. A basic aspect of this problem is finding a criterion for determining when 
a matrix has a pair of pure imaginary eigenvalues. A complete algebraic solution to this 
problem is more complex than finding a criterion for zero eigenvalues. The set of matrices 
with pure imaginary eigenvalues does not define an algebraic hypersurface in the space of 
matrices, but rather a semialgebraic set. Let us see what this means in concrete terms 
for 2x2 matrices. A 2 x 2 matrix has pure imaginary eigenvalues if and only if its trace 
vanishes and its determinant is positive. Thus the set of 2 x 2 matrices with pure imaginary 
eigenvalues is specified by one equation and one inequality on the coefficients of the matrix. 
If we drop the inequality, then we find matrices with real eigenvalues of equal magnitude 
and opposite sign as well as matrices with pure imaginary eigenvalues. 

There are classical algebraic theories that address the question of when a matrix has a 
pair of pure imaginary eigenvalues. These theories are now best known in the form of the 
Routh-Hurwitz criteria that a matrix have all of its eigenvalues in the left half plane. We 
have applied these theories and continuation strategies in order to implement algorithms 
for computing Hopf bifurcations. Here we give a brief indication of the algebraic foundation 
for our algorithms that give criteria for a square matrix to have a pair of pure imaginary 
eigenvalues. For the remainder of this section, we let A denote an n x n real matrix and 
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P(A) = det{\I — A) denote its characteristic polynomial. The roots of P are the eigenvalues 
of A. There are two approaches to our algebraic problem: To give criteria for the polynomial 
P to have a pair of pure imaginary roots, or, alternatively, to give criteria for A to have a 
pair of pure imaginary eigenvalues by performing algebraic transformations directly on A. 

We make the basic observation that if a real polynomial P has a pure imaginary root A, 
then — A is also a root. Thus a necessary condition for P to have a pair of pure imaginary 
roots is that the polynomials P and Q{x) = P{—x) have a common root. The Sylvester 
resultant is a function of the coefficients of P and Q that vanishes if and only if P and Q 
share a common root. Applied in the context of our problem, we first reduce the size of 
the problem by observing that the polynomials R = P + Q and S = P — Q are even and 
odd, respectively. Moreover, P and Q have a common root if and only if R and S have a 
common root. We can write R{x) = R{x'^) and S{x) = xS{x'^) with R and S polynomials 
of degree approximately n/2. Therefore, the Sylvester resultant of R and S gives a function 
that vanishes if and only if P and Q have a common root or, equivalently, P has two roots 
whose sum is zero. If we denote 

P(A) = Co + CjA + • • • + c„_iA""' + A" , 



and n is even, the Sylvester resultant is the determinant of the matrix 



Co C2 • • • 

Co C2 
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V Ci C3 ••• 

while if n is odd, it is the determinant of the matrix 
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The resultant encodes the outcome of applying the Euclidean algorithm to a pair of 
polynomials. If two polynomials have a common root, then the Euclidean algorithm yields 

the greatest common denominator of the two polynomials. The coefficients of this GCD 
can be expressed as determinants in the coefficients of the two polynomials. These determi- 
nants are the subresultants of the two polynomials. These resultants give explicit functions 
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n 


P(A) = Co + ciA + . . . + c„_iA"-i + A" 


2 


ci = , Co > 


Q 
o 


CQ — ClC2 — U , Cl ^ u 


4 


coc| - C1C2C3 + cf = , C1C3 > 


5 


(C2 - C3C4)(ciC2 - C0C3) + CiC4(ciC4 - 2co) + cg = , 
(C2 - C3C4)(co - C1C4) > 


6 


Coci(coC5 - C2C3) + Cic|(c| - C0C4) + Ci(cf + C0C3C5) + CiC5(coC3 - 2ciC2) + 
(C4C5 - C3)(coc| - C0C1C5 + cfc4 - C1C2C3) = , 
(C1C3 + Coci - CiC4C5)(c§ - C1C5 + C2C| - C3C4C5) > 



Table 1: Conditions for P{\) to have a pair of pure imaginary roots. 



that define the locus of matrices that have eigenvalues whose sum is zero. They do not 
determine whether the roots are real or complex. For this purpose, one can use the theory 
of subresultants |l49[| . We observe that a common root of P and Q is imaginary if and only 
if the corresponding common root of R and 5 is negative. Therefore, we can give explicit 
inequalities in terms of the subresultants that determine whether a polynomial with a single, 
simple pair of roots whose sum is zero has a pair of pure imaginary roots. Table |l| gives a 
list of the equations and inequalities that determine the polynomials P of degree at most 
six with these properties. 

The computation of characteristic polynomials of square matrices is computationally 
expensive and prone to numerical errors for some types of matrices. Thus, we seek methods 
that allow us to determine whether a matrix has a pair of pure imaginary eigenvalues 
without computing the characteristic polynomial as an intermediate step. This can be done 
through the definition of appropriate Kronecker or tensor products. The basic algebraic 
idea is that there are transformations of matrices that produce matrices whose eigenvalues 
are functions of the eigenvalues of the original matrix. For example, if A and B are square 
matrices, then the eigenvalues of the tensor product A<^ B are the pairwise products of the 
eigenvalues of A and those of B. If we form the transformation T = I<^A + AfSiI with 
/ the n X n identity matrix, then the eigenvalues of T are sums of pairs of eigenvalues of 
A. To remove the redundancy associated with having both pairs Aj + Aj and Xj + Aj as 
eigenvalues, we use a skew-symmetric version construction. The bialternate product AQ I 
of A is the n(n — l)/2 x n(n — l)/2 matrix defined by 



(2A0 1„; 



{pq,rs} 



-iA)p,s ifr = g 

{A)p^r if r 7^ p and s = 

(^)p,p + {A)q^g if r = p and s = 

{A)q^s if r = p and s / 

-iA)q^r iis=p 

otherwise 



where the rows are labeled lexicographically by for {p = 2, . . . ,n;q = 1, ... ,p — 1) and 
the columns likewise by rs for (r = 2, . . . , n; s = 1, . . . , r — 1) [Q. The eigenvalues of 
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A Q I are pairwise sums Aj + Xj with i < j of the eigenvalues of A. Thus, a necessary 
condition for A to have a pair of pure imaginary eigenvalues is that ^ / be singular. 
This singularity can be tested by computing the determinant oi A Q I, but we can also 
use other algorithms (such as singular value decomposition) from numerical linear algebra 
that give more accurate and robust tests for the singularity of A. In our ongoing work, 
we are implementing algorithms for tracking Hopf bifurcations based upon these algebraic 
methods and standard continuation methods. 

4 Homoclinic Bifurcations and Their Computation 

Bifurcation theory describes qualitative changes in phase portraits that occur as parameters 
are varied in the definition of a dynamical system. For dynamical systems defined by vector 
fields on i?", one has a system of equations of the form x = fx{x) with A S denoting 
a /c-dimensional vector of parameters. A rough classification of bifurcations distinguishes 
between local and global bifurcations, but it is difficult to make this distinction precise. 
Heuristically, the dynamics of local bifurcations are determined by information contained 
in the germ of the Taylor series of / or a Poincare return map at a point, whereas the 
dynamics of global bifurcations require information about the vector field along an entire 
(non-periodic) trajectory. While the theory of bifurcations in multi-parameter families is 
far from complete, the theory of global bifurcations is more fragmentary than that of local 
bifurcations. A number of "codimension 2" global bifurcations have been studied, but 
there has not been an attempt to construct a synthesis of these studies that portrays a 
systematic view of the different cases and phenomena that occur. The lectures that were 
given in Montreal at this NATO sponsored summer school touched upon these matters, 
but these notes go much farther towards the construction of such a synthesis. We focus our 
attention on global bifurcations that involve homoclinic or heteroclinic orbits of equilibrium 
points of a vector field. We shall call homoclinic and heteroclinic orbits connecting orbits 
in order to have a common name for the two. 

Silnikov seems to be the principal originator of the strategy we adopt to study global 
bifurcations with connecting orbits in higher dimensional vector fields. For planar vector 
fields, the techniques are older and asymptotic methods for studying global bifurcations 
in planar vector fields were well developed by the 1920's (Dulac [^]). The fundamental 
idea is that the recurrent behavior near a connecting orbit should be studied in a fashion 
similar to that used in studying periodic orbits via a Poincare return map. In particular, 
codimension one cross sections to the flow are introduced, and the return map of these 
cross sections is studied. There are some additional complications in the study of connect- 
ing orbits compared to that of periodic orbits which significantly complicate the analysis. 
We mention three of these. First, the discrete maps that describe fiow past an equilibrium 
point are singular. In the simplest cases of fiow past non-resonant equilibrium points, these 
"passage" or "correspondence" maps have singularities of the form where a is not an 
integer and may be complex if the equilibrium has complex eigenvalues. These singularities 
lead to analysis that is much more complicated and intricate than that associated with the 
return map of a periodic orbit. As the work presented by Ilyashenko and Ecalle at this 
summer school demonstrated, this analysis is frightfully complicated even for planar vector 
fields. The second complication arising from connecting orbits is that there are discontinu- 
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ities in the return maps of cross sections that are associated with connecting orbits. The 
bifurcations and attractors that appear in the Lorenz system [50| give a vivid example of 
the consequences of these discontinuities and have been described by Guckenheimer and 
Wilhams |37|. These discontinuities force one to look at dynamical systems that are built 
from multiple pieces rather than studying the iterates of a single continuous mapping. The 
third complication is that the list of cases of qualitatively different type is substantially 
longer than with local bifurcations. Tools for treating all of the cases in a single analy- 
sis are lacking, so the construction of a comprehensive and complete theory seems like a 
daunting task. 

These lectures adopt a superficial view of the mathematical technicalities associated 
with global bifurcations. In constructing return maps for a flow along connecting orbits, one 
would like the simplest possible analytical expressions for these return maps. Normal form 
theory explores how coordinate changes can be used to simplify these analytical expressions, 
but the theory produces a large and intricate story whose details can mask many of the 
phenomena that we want to study. The idea of Silnikov was to relax the requirement 
that an exact normal form be used to describe passage past an equilibrium and to use an 
approximation to the normal form instead. Asymptotic expansions for the flows and passage 
maps can be constructed, and one can begin the analysis of bifurcations with connecting 
orbits by using the leading term, or terms, of these asymptotic expansions. A similar 
approach to the global part of the flow producing returns to a cross section can also be used, 
but the asymptotic expansions of these smooth maps are given simply by Taylor series. A 
matching procedure can then be used to represent the full return map as a composition 
of singular and regular maps that come from passage past equilibrium and parts of the 
flow that are non-singular. In constructing maps with such matching procedures, we must 
remember that there may be discontinuities that lead to the examination of several different 
sequences of compositions. Following an analysis of return maps built from approximations, 
we can seek to determine the structural stability of the systems. It is unreasonable to expect 
that each type of bifurcation will have a structurally stable unfolding. We are plunged into 
the morass of considerations that result from the fact that structurally stable vector fields 
are not dense in the space of vector fields. The most that we try to do is to explore in each 
case which dynamical features of the bifurcation do persist under perturbations. Note that 
even specifying which perturbations are to be allowed in a theory based upon the return 
maps of connecting orbits is a tricky matter due to the singularities of the maps. Rather 
than engaging the reader in an extensive discussion of these problems, we proceed past them 
while erecting signposts that point in the direction of unresolved and incomplete technical 
matters. 



4.1 Generic Homoclinic Orbits 

Let F : i?" BP' be a sufficiently smooth, vector field with a homoclinic orbit to an 
equilibrium point at 0. We shall denote the flow of F by (pt and the homoclinic orbit by 
7(t). Let 

r = {7(t) : t G i?} where lim -f{t) = and 7(t) / 0. 

We give a list of conditions that generic systems with a homoclinic orbit satisfy. These 
conditions help ensure that bifurcations from the homoclinic orbit have as simple a structure 
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as possible. In classifying codimension two bifurcations of homoclinic orbits, we can then 
look at events that cause the failure of one of these genericity conditions and examine the 
effect on bifurcations in two parameter families. 

The first condition is that the equilibrium be hyperbolic and non-resonant. The non- 
resonance conditions that are important are primarily ones that involve the stable and 
unstable eigenvalues of smallest magnitude. We call these eigenvalues the principal eigen- 
values of the equilibrium. More precisely, assume the linearization DxF(0) has k eigenvalues 
in the left half plane and n — A; in the right half plane. Assume there exists a real principal 
eigenvalue in the stable manifold denoted or a complex conjugate pair Xs,Xs with the 
property that if A is another stable eigenvalue of Dj;F(0), then Re(A) < Rc(As) < 0. Simi- 
larly, assume there exists a real principal eigenvalue in the unstable manifold denoted A^ or 
a complex conjugate pair A^, Xu with the property that if A is another unstable eigenvalue 
of D^F{0), then < Re(A„) < Re(A). 

At points in the stable and unstable manifolds of 0, there are filtrations of the tan- 
gent bundle associated with the exponential growth or decay of vectors as the trajectory 
approaches the origin. For a point p of the stable manifold, we are interested in three 
subspaces that we denote E^~^{p) D E^{p) D E^^{p) and call the stable plus weak unstable 
manifold, the stable manifold and the strong stable manifold, respectively. These are defined 
via the variational equation for the trajectory through the point p given by 

Let Vt{v;p) denote the solution to this linear equation with initial condition Vo{v;p) = v e 
TpR^. Let Vs and be numbers such that 

Re(A) <Us < Re(As) < or < Re(Au) < f « < Re(A). 

for any eigenvalue A of DxF{0) which is not a principal eigenvalue. Then 

E'+{p) = {ve TpR^l lim^^+oo e-''-*\Vtiv-p)\ = 0} 

E'{p) = {ve TpR^\ lim^^+oo \Vt{v;p)\ = 0} 
E''{p) = {v& TpR^l limt^+oo e-'^'\Vt{v;p)\ = 0} . 

There are analogous definitions of the unstable plus weak stable manifold, the unstable 
manifold and the strong unstable manifold of a point p in the unstable manifold of 0: 

E'^+ip) = {ve Tpi?"! lim^^.oo e-'''\Vt{v;p)\ = 0} 
= {^J G Tpi?"! limt^_oo \Vt{v;p)\ = 0} 
= {ve TpR^'l limt^_oo e-^"*|Vt(i;;p)| = 0} . 

Note that there is no apparent relationship between these sets of manifolds at a point p of 
a homoclinic orbit other than that the vector field F{p) belongs to both E^(p) and E'^{p). 

Our next requirement for a generic homoclinic orbit involves the direction of approach 
of a homoclinic orbit to as t — ±00. Associated with the origin itself is a strong stable 

manifold W^^ consisting of points p G W^, the stable manifold of the origin as defined by 
the Stable Manifold Theorem, for which the vector field F(p) lies in E^^{p). Note that 
TqW = E''{0) just as ToW = E'{0). Now W is a proper submanifold of W of 
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codimension one or two depending upon whether the principal stable eigenvalue is real or 
complex. In either case, almost all trajectories in lie in Similarly, almost all 

trajectories in lie in W'^\W'^^. Our next requirement for a generic homoclinic orbit is 
that a point p on F satisfy 

In words, (H3) is the statement that the homoclinic orbit approach the origin in the direc- 
tions of the principal eigenvectors. This may alternatively be expressed by requiring that a 
point p on r satisfy 

F{p) E'', F{p) . 

The final condition that we impose upon generic homoclinic orbits involves the derivative 
of the flow along the homoclinic orbit. We state the condition in terms of the stable plus 
weak unstable and the unstable plus weak stable manifolds. Depending upon the types of 
the principal eigenvalues, the sum of the dimensions of E^^(p) and E^~^{p) is n + 2, n + 3, 
or n + 4. We require that the intersection E^^{p) n E^~^{p) be transverse and, furthermore, 
that E'+{p) n E'^'^{p) = E^'ip) n E^'+ip) = {0}. 

Let us examine the meaning of this final condition in the case of real eigenvalues. Assume 
that As, \u G R- Then dim(£'*'*) = m — 1 and dim(i?"") = n — 1. Along T, the intersection 
W^~^ n W^~^ is a two-dimensional bundle that contains the vector field. At the origin, 
this bundle approaches the plane spanned by the principal eigenvectors, both for t — > 
-1-cxD and for t —oo. We can think of this bundle as a "ribbon" along the homoclinic 
orbit that defines the behavior of the system in the "weak" directions that determine the 
primary structure of the orbit. Taking the closure at the origin of the bundle along the 
homoclinic orbit, we obtain a bundle of planes along a simple closed curve. This bundle 
is either orientable or non-orientable. We distinguish these cases by calling them twisted 
and untwisted homoclinic orbits. Twisted homoclinic orbits cannot occur for vector fields on 
orientable two-dimensional manifolds. Both twist types of homoclinic orbits are represented 
pictorially in Figure ^. 

We will call a homoclinic orbit with real principal eigenvalues a binodal homoclinic orbit. 
For the generic binodal homoclinic orbit, there are no additional interesting dynamical 
structures in a sufficiently small neighborhood of the homoclinic orbit ||6^. This is not 
necessarily true of generic homoclinic orbits with a complex principal eigenvalue. When 
exactly one of the principal eigenvalues is complex, we will call the homoclinic orbit a 



unifocal homoclinic orbit, and the generic case has been studied by Silnikov pO| , 62, 63] and 



others [24, 25, ESl, p&. The following two theorems describe the dynamical structures close 



to a generic unifocal homoclinic orbit, assuming that Xu ^ R and E C. 

Theorem 4.1 // |A„/i?e(As)| > 1 (Silnikov condition) then there exist horseshoes in every 
neighborhood ofT. 



Theorem 4.2 // \Xu/ Re{Xs)\ < 1 then there is a neighborhood of F which contains no 
periodic orbits. 

The non-resonance condition on the principal eigenvalues is given by |At(/Re(As)| ^ 1. 
We may look at further conditions on the eigenvalues and determine a second genericity 
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requirement: |Au/Re(As)| 7^ 2. If 1 < |AM/Re(As)| < 2, then the hnear flow at the origin is 
contracting and the horseshoes are attracting, while for |A„/Re(As)| > 2 the linear flow at 
the origin is expanding. 

Finally we consider the generic homoclinic orbit with complex principal eigenvalues 
which is called a bifocal homoclinic orbit. Figure § depicts the unifocal and bifocal homo- 



clinic orbits. The bifocal homoclinic orbit has been studied by Silnikov |61, 33 1, Glendinning 



|28|, and Fowler and Sparrow |pO[. They prove the following theorem: 



Theorem 4.3 There exist horseshoes in every neighborhood of a generic bifocal homoclinic 
orbit. 



4.2 Silnikov Coordinates and Local Normal Forms 

The unfoldings of the bifurcations involving connecting orbits can be very complicated. As 
we stated above, certain types of generic homoclinic orbits are embedded in much more 
complex dynamical structures. A full description of the unfolding of these orbits entails a 
comprehensive analysis of the horseshoes that are created or destroyed as a parameter is 
varied. These unfoldings have been studied for unifocal homoclinic orbits [p4|, p5|, [29|] and 



bifocal ones |28[]. We do not undertake such a monumentally complicated task, but we 
do want to describe maps that give an approximation to the return maps associated with 
homoclinic orbits. The general procedure we employ for doing this involves a decomposition 
of the return map for a homoclinic orbit into two parts: one describing the "local" part 
of the flow past the equilibrium point in the homoclinic orbit and the other describing 
the "global" portion of the flow outside of this neighborhood. In cases of singular cycles 
containing more than one connecting orbit, we make a decomposition into more pieces, but 
the principle is the same. In the decompositions that we use, we try to choose coordinates 
in a manner that simplifies the analytical expressions of the vector fields. 

The simplest flows near an equilibrium point are linear. The question as to whether 
coordinate transformations near an equilibrium point can be found that linearize a vector 
field in a neighborhood of the equilibrium has been studied systematically for over a century, 
starting with Poincare's dissertation. The answers to the question are complicated. See 
Arnold for an extensive summary of what is known about the linearization problem. 
Here we shall use a few bits of this theory. Resonance conditions on the eigenvalues provide 
the most elementary obstructions to linearization. A resonance condition of order k is 
expressed in terms of eigenvalues Xj hy Xi = J2 o-j^j with the Oj nonnegative integers whose 
sum is A;. A vector field with an equilibrium point that satisfies a resonance condition of 
order k usually cannot be linearized by a coordinate transformation. Nonetheless, there 
are resonant normal forms for the equilibrium with the property that there are smooth 
coordinate changes that transform the degree I Taylor series of the vector field at the 
equilibrium to its resonant normal form. We can hope in these situations that it is feasible 
to describe explicitly the flow of the resonant normal forms truncated to degree I and that 
these flows serve as good approximations to the flow of the original vector field in the vicinity 
of the equilibrium point. Here the classical theory breaks down and does not provide a good 
solution to the question as to when the passage map of two flows near an equilibrium are 
good approximations to one another. 
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Consider a linear vector field x = f[x) with a hyperbolic equilibrium point at the origin 
with k stable eigenvalues and n — k unstable eigenvalues. Choose a coordinate system so 
that the stable and unstable manifolds and of the origin lie in coordinate planes. Let 
Uq be a neighborhood of the origin that is the product of balls and of radius r in the 
stable and unstable manifolds and W^. The boundary of x B^ is dB^ x B^yjB^. x dB"^. 
Since the vector field / is linear, the motion of points is the superposition of motions along 
the stable and unstable manifolds. Furthermore, the distance to the origin of points in the 
stable manifold decreases monotonically, while the distance to the origin of points in the 
unstable manifold increases monotonically, in a well chosen coordinate system. Therefore 
the passage map of Uq will be defined as a map (f) : dB^. x 5," B'^. x 55". With this choice 
of neighborhood of the origin, it is difficult to determine an explicit expression for the time 
of flight for a trajectory to reach the outgoing boundary of Uq and hence to obtain a formula 
for (f). On the other hand, if we choose different neighborhoods of the origin, then we can 
sometimes compute the exit time from the neighborhood. In the case of real eigenvalues, 
we find that there are power law singularities; complex eigenvalues produce singularities 
with other elementary functions. Explicit examples are computed later. 

For some problems of higher codimension, we shall need to compute the passage maps 
of equilibrium points that are either resonant or nonhyperbolic. In these cases, we shall 
still seek coordinate systems and approximations for which there arc explicit integrals for 
the local normal forms at the equilibrium point. Cases with this property are the only 
ones considered in this paper, though there are higher codimension problems for which 
the normal forms are not integrable. In finding passage maps for these resonant cases, we 
would like to obtain formulas that remain valid when we perturb the equilibrium to make 
it generic. This requires some additional care beyond merely solving the explicit flows for 
the passage map at the resonant equilibrium. 

The Silnikov procedure is to combine the local passage maps near equilibria with nonsin- 
gular maps that describe the flow between cross sections around the portions of connecting 
orbits that do not contain equilibria. In the case of a homoclinic orbit, these cross sec- 
tions can be taken to be portions of the boundary of the neighborhood Uq that was used to 
construct the passage map past the equilibrium. As with the passage maps, we seek approx- 
imations for these "global" maps between cross sections. Since the maps are nonsingular, 
they can be approximated by truncating their Taylor series. Generally, we start with affine 
approximations to the transformations. If these are inadequate to obtain the (structural) 
stability results we seek, then higher degree approximations are used. The return map for 
a cross section can be obtained by composing these global maps with the local passage 
maps. In carrying through this composition, care must be taken with understanding the 
domains on which the maps are defined. When there are multiple connecting orbits, there 
is the additional possibility that flow in different directions away from an equilibrium may 
produce structures that require following different patterns of return to a neighborhood of 
the equilibrium. 

4.3 Bifurcations From a Homoclinic Orbit: What Can We Study? 

The difficulties of proving rigorous theorems about the unfoldings of bifurcations are formid- 
able. Nonetheless, we would like to prove as much about phenomena that are consequences 
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of bifurcations of homoclinic and heteroclinic orbits as we can. Indeed, we would like to 
capture the major dynamical events that occur in the unfolding despite the fact that struc- 
tural stability of the unfoldings will often fail. To make the theory "local" to the bifurcating 
orbits, we restrict attention to an arbitrarily small neighborhood of these orbits. Recurrent 
behavior that involves trajectories that leave such a neighborhood will be discussed sepa- 
rately. Here we describe three types of structures that can occur in small neighborhoods of 
bifurcating connecting orbits. 

The first type of dynamical structure that can be involved at a homoclinic bifurcation 
is a periodic orbit. Generic bifurcations of binodal homoclinic orbits produce periodic 
orbits whose stability is determined by the sum of the principal eigenvalues at the saddle. 
In higher codimensions, multiple periodic orbits can bifurcate. In some cases (notably, the 
"gluing bifurcation" 1^^), these orbits and their patterns of bifurcations can be complicated. 
The second type of dynamical structure produced in bifurcations of connecting orbits is 
typified by the homoclinic bifurcation in the Lorenz system ||5^. In this system, for reasons 
of symmetry, a pair of homoclinic orbits are formed. As they bifurcate, a horseshoe is 
immediately created. The same phenomenon occurs in higher codimension bifurcations 
without symmetry. Third, generic homoclinic orbits with complex values may be adjacent 
to horseshoes. As described by Silnikov [^], homoclinic orbits to equilibria satisfying 
appropriate conditions on their eigenvalues occur only in the closure of horseshoes. As the 
homoclinic orbit is unfolded, the horseshoes closest to the homoclinic orbit can be destroyed. 

The singularity theory for mappings between two spaces is much "cleaner" than the 
bifurcation theory of flows (or iterated mappings). There are a number of reasons for these 
differences, one being the complex dynamical structures that occur in unfoldings of certain 
bifurcations. If there are horseshoes that are created or destroyed as part of a bifurcation, 
then the unfolding of these bifurcations will involve all of the associated complications. 
These have been studied most thoroughly in the context of the Henon mapping 52 1. 
There are an infinite number of bifurcations of periodic orbits that are part of the creation 
of horseshoes. If there is a single unstable eigenvalue and volume contraction in these flows, 
there are also phenomena such as infinitely many (Newhouse) sinks and the occurrence 
of nonhyperbolic strange attractors |3^. We do not want to become enmeshed in these 
details. To a large extent they appear to be subsidiary to the fact that horseshoes are 
created, and we expect to add little to the general discussion of the processes associated 
with the creation and destruction of the horseshoes. Our focus in dealing with horseshoes 
will be to demonstrate that a return map in an unfolding satisfies the conditions that 
guarantee the existence of horseshoes. We will not generally explore whether the horseshoe 
is part of a larger invariant set or investigate the presence of nonhyperbolic invariant sets 
in the unfoldings. 



4.4 Codimension Two Bifurcations of Connecting Orbits 

There are many types of codimension two bifurcations of connecting orbits. Failure of 
one of the conditions that characterize a generic homoclinic orbit will lead to a degenerate 
bifurcation. These can be classified into four groups: 

1. Eigenvalue degeneracies 
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2. Degenerate approach 

3. Degenerate twist 

4. Multiple connecting orbits 

Each of these types can be further subdivided. For example, the eigenvalue degeneracies 
may be due to a single zero eigenvalue, a pair of pure imaginary eigenvalues, equal magnitude 
of the real parts of the principal eigenvalues, or equal magnitude of a real principal eigenvalue 
with the sum of a pair of imaginary principal eigenvalues. Furthermore, one has a different 
analysis depending on whether the homoclinic orbit is nodal, unifocal or bifocal, and, in 
the nodal case, the twist associated with a resonance. In all of the cases that have been 
studied, the introduction of Silnikov coordinates and the study of return maps built from 
these coordinates is a major portion of the analysis of the bifurcation. Rigorous results that 
go beyond the description of model systems tend to be very difficult to formulate and prove. 
When there are horseshoes associated with these bifurcations, even the phenomenology 
associated with the Silnikov approximations tends to be incomplete. 

Below we discuss a bifurcation with an eigenvalue degeneracy, one with a degenerate 
approach, and another with multiple connecting orbits, describing each in terms of Silnikov 
approximations. Additionally, within each category, we give a survey of other work known 
to us. For the case of degenerate approach, we refer the reader to the work of Terman |65|. 



4.5 Eigenvalue Degeneracies 

There are a number of different eigenvalue degeneracies which may occur. Non-hyperbolic 
equilibrium point degeneracies have been studied by, in the nodal case, Deng Luk'yanov 
| |51| and Schecter ||5^; in the case of unifocal homoclinic orbits with a nonhyperbolic saddle, 
Belyakov j^; and, in the case of unifocal homoclinic orbits with a nonhyperbolic focus, 
Argoul, Arneodo and Richetti |l|, |5^, Bosch and Simo ||8|, Gaspard and Wang [^], and 
Hirschberg and Knobloch ||4^. Here we will look at vector fields that have a homoclinic 
orbit with real principal eigenvalues of equal magnitude. The homoclinic orbit may be 
twisted or untwisted, and the bifurcation is different, depending on the case. This has been 



studied by a number of different authors, including Glendinning |27], Kokubu |47], and 
Chow, Deng, and Fiedler @. 

A normal form for a planar saddle with a 1:1 resonance is given by the system 

X = — x(l + xy) 

y = y 

An unfolding of this system allows the ratio of eigenvalues at the origin to vary: 

X = —x{l + A + xy) 

y = y 



To use Silnikov coordinates for these systems we seek to integrate the normal form explicitly 
and solve for the passage map past the origin. The trajectory with initial conditions {xo,yo) 
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e-(i+^)« 



yit) = yoe* 

and the passage map from the hne x = 1 to the hne y = 1 has a time of flight t = — ln(y) 
and an expression 

^ = V'(y) = mo^^ • 

When A — > 0, this map becomes 

y 



X = i}{y) 



l-ylny 



The leading order term in Equation Q is and is adequate for the calculation of the 

unfolding of the codimension two bifurcations of the homoclinic orbit. The global portion 
of the Silnikov approximation to the return map is an affine transformation y = ax + fi. 
We assume that a 7^ ±1 and that /i is the second parameter of the unfolding. Thus the 
approximate return map is </>(y) = ay^^^^^ + fi. 

If the homoclinic orbit is untwisted, then a > 0. If the homoclinic orbit is twisted, then 
o < 0. The domain of (/> is a segment with endpoint at y = 0. We analyze (p first in the 
untwisted case. There is a periodic orbit close to the homoclinic orbit if (j) has a fixed point. 
The periodic orbit is degenerate if its fixed point p satisfies (l)'ip) = 1- This happens when 
and a{l + X)p^ = 1. If we eliminate p from this pair of equations, we obtain 

^ = A(l + A)-(i+i/^)a-V^ 

As A — > 0, this gives /x ~ Xe^^a^^^^ . In order that // be small, we need to choose the sign 
of A so that A In a > 0. Thus A is positive if a > 1 and A is negative if < a < 1. Observe 
that the curve in the (A, p) plane along which nonhyperbolic periodic orbits occur has a 
flat tangency with the curve /x = along which there are homoclinic orbits. In the thin 
wedge between the two, there are two periodic orbits near the location of the degenerate 
homoclinic orbit. See Figure ^. 

The analysis of the twisted case is similar, but there is indeed a new twist. The return 
map has the form described above, but a < 0. This makes the return map a decreasing 
function of y. As a result, there are two kinds of periodic orbits and two kinds of homoclinic 
orbits in the unfolding: "once rounding" and "twice rounding". The once rounding loops 
are twisted while the twice rounding loops are untwisted. The once rounding periodic 
orbits have a negative principal characteristic multipler while the twice rounding orbits 
have a positive multiplier. The transition between the two is given by a period doubling 
bifurcation. Thus there are three different types of bifurcations that occur in the unfolding 
of the twisted resonant loop: once rounding homoclinic orbits, twice rounding homoclinic 
orbits and period doubling bifurcations of periodic orbits. These can be computed in terms 
of the return map once rounding homoclinic orbits are given 

by = 0. The twice rounding homoclinic orbits are given by parameters for which = 



Dynamical Systems: Some Computational Problems 



17 



0^(0) = (f>{fi) = afi^^^^^ + fx. Solving this equation for /x yields fi = |a|~^/'*'. The period 
doubling bifurcations occur at fixed points of cj) for which the return map has derivative —1: 
a(l + X)p^ = — 1 and /i = p — ap^^'^^\ This yields 

As A — > 0, observe that the ratio of the values of /j, that produce period doubling and twice 
rounding homoclinic orbits approaches 2/e. This ratio is a "universal" property of this 
codimension two bifurcation. As with the untwisted resonant bifurcation, the bifurcation 
curves have a flat tangency. See Figure ^. 

4.6 Inclination Degeneracies 

When the principal eigenvalues of a homoclinic orbit are real, but the homoclinic orbit is 
neither twisted nor untwisted, then there is an inclination degeneracy. This situation has 
been studied by Deng [|l^, and there are many subcases that lead to different bifurcation 
diagrams. Here we describe one of these cases, illustrating that horseshoes can occur in 
the unfolding of this codimension two global bifurcation. The example also illustrates a 
circumstance in which an approximate return map that comes from the composition of 
an affine map with a local passage transformation does not capture all of the important 
dynamical behavior. 

Three is the lowest dimension in which one can construct a homoclinic orbit with de- 
generate twist. We consider a two parameter family of vector fields that are linear in a 
neighborhood of the origin, with two real stable eigenvalues Ay < A^; < and an unstable 
real eigenvalue A^ > 0. We assume that there is a homoclinic orbit that approaches the 
origin along the x axis in forward time and along the z axis in backward time. We shall 
assume further that the return map for the orbit is degenerate in the twisting of the normal 
bundle around the homoclinic orbit. This means that a vector, pointing in the direction 
of the principal stable eigenvector as the homoclinic orbit leaves the origin, returns in the 
direction of the strong stable manifold (rather than the generic behavior of returning in the 
direction of the unstable manifold). 

Let us view the situation in terms of a Silnikov approximation. The passage map near 
the origin from the cross section S defined by x = 1 to the cross section z = 1 will be given 
by 




where 

a = -Xx/>^z, /? = -Ay/Az, < a < (3 . 

Here and below, if z < 0, we denote z°^ = — (— The global return from ^ = 1 to x = 1 
will be approximated by 

( y \ _ / au + bv \ 

[z )-^[v )-[x )^[ -^,u + u^ + v ) ■ 
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Then the return map (f) is approximated by 

_ / c + az° + byz^ \ 
~ \ A - ;uz" + + yzl^ ) ■ 

Here the parameter A determines whether a homochnic connection occurs, while the sign of 
/i determines whether the connection is twisted or untwisted. Note that the approximation 
has included a quadratic term in 7. Without this term, the plane spanned by the principal 
eigenvectors would return entirely within the stable manifold, a highly degenerate situation. 

There are now several cases determined by the relative magnitude of the exponents a 
and (3. The case we shall examine is the one in which we assume that 2a < (3 and a < 1/2. 
The first of these assumptions implies that \yz^\ = o{z'^'^) and the second leads to complex 
dynamical behavior involving the stretching in the z direction. We shall not try to give a 
complete analysis of even the specific family of maps defined by (f), but we shall show that 
there are horseshoes that appear in some regions of the parameter space for this family. 

The assumption 2a < (3 leads to an approximation of by a one-dimensional mapping. 
The nature of the approximation will be examined below. Consider the one-dimensional 
mapping 0(z) = A — fiz" + z^" obtained by ignoring the y coordinate in the definition of (p. 
This is a unimodal map, and we look for ranges of parameter values for which this mapping 
will have a hyperbolic invariant set. If A = S/x^/lG and /x > 0, = {z" —ij,/4){z°' — 3^/4). 
Let 

A = 0, B = (;u/4)^/", C = (3;u/4)i/", D = /i^/" , 

so A < B < C < D. Then (j){A) = ^{D) = 2,iJ? /IQ > D since a < 1/2, and ^{B) = 
0(C7) = A. Now W{B)\ = W{C)\ = 2a(-/x/4)2-V« > 1 for smah. Since ^{[A,B]) D 
[A,D], 4){[C,D]) D [A,D], and (p > I on [A, B] and [C, D], we conclude that the map (p has 
a hyperbolic invariant set in the interval [A,D) = [0, /i^/"). 

In the cross section E, consider the rectangle = [c—2^^/",c+2yL^/°'\ x [—2fi^/", 2/x^/"]. 
If we set A = S/z^/lG and let ^ — > from below, then the return maps of the rectangles 

can be rescaled so that they tend to a map of rank 1. The rescaled maps have a 
vertical coordinate that behaves as a rescaled version of (p. This situation is analogous to 
the behavior of the Henon map |4C] in the limit as the map tends to a map of rank one. The 



hyperbolic invariant set will persist as a horseshoe as we leave the singular limit /x — 0. The 
analysis of Benedicks and Carleson |p, extended by Mora and Viana [^] indicates further 
that as A varies, we will expect to encounter chaotic attractors for the return map (p. 



4.7 Multiple Connections 

Consider a vector field in a generic two-parameter family that has an equilibrium point with 
a single unstable eigenvalue, real non-resonant principal eigenvalues and a pair of generic 
homoclinic orbits formed from the two unstable separatrices of the saddle point. We would 
like to analyze the unfoldings of such systems. In the planar case, this determines the 
geometry of the unfolding of the family and is described below. On the other hand, there are 
several choices that occur if the vector field has dimension at least three. In particular, if the 
two homoclinic orbits approach the equilibrium from the same direction along the principal 
stable eigenvector, then there are cases that lead to either geometric Lorenz attractors 
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or to complex periodic orbits formed by the gluing bifurcations described by Gambaudo, 
Glendinning and Tresser [^]. We recall bits of their analysis below. 

Let X be a planar vector field with a non-resonant hyperbolic equilibrium point at the 
origin and a pair of orbits homoclinic to the origin. Assume that X is contained in a generic 
two parameter family of vector fields. There are two topologically distinct configurations 
depending upon whether one homoclinic orbit is contained in the other, but there is little 
difference in the analysis of the two cases. Pick a pair of cross sections Si,E2 to the two 
stable separatrices and consider the return map </> : S1UE2 — > S1US2. This return map will 
have discontinuities at the intersections of Si and S2 with the stable manifold of 0. Thus we 
can think of (j) as constructed from four maps (pij : Sj T,j . These maps fit together so that 
the ranges of cpij and (j)2j are contiguous and do not overlap. At the common endpoint of 
the images, there is a singularity of the form with a 7^ 1. By reversing time if necessary, 
we may assume that the magnitude of the stable eigenvalue of is larger than the unstable 
eigenvalue, implying a > 1 and that the derivatives of the (pij are smaller than 1 near the 
stable manifold of 0. Furthermore, we may choose orientations of Si and S2 so that the 
maps (pij are increasing. Then each (pij or iterates of the (pij in a connected domain can have 
at most one fixed point, and that point will be stable. Determining which types of periodic 
orbits bifurcate from the double loop becomes a matter of determining which configurations 
of fixed points for iterates of the maps (pij can occur. The cross sections and return maps 
are pictured in Figure |^. 

Denote once again = — (— x)" for x < and choose coordinates on cross sections 
that are centered on the stable manifolds. The Silnikov approximation for the return maps 
have the form (pij{x) = ajx" + 6j, but one must remember that there are two components 
to the domain of the return map and (j)i2 and (j)2i map one domain to the other. Here 61 
and 62 can be regarded as the unfolding parameters of the system. 

It is important here that the maps and (j)2j have contiguous images. There are three 
kinds of periodic orbits for the return maps that correspond to periodic orbits of the flow: 
fixed points of 0^ and fixed points of a composition (^12 o (f)2i or (f)2i o (j)i2. The fixed points 
of correspond to periodic orbits that lie close to one of the homoclinic orbits, while the 
fixed points of 0i2 o (j)2i and </>2i o ^12 correspond to periodic orbits that lie close to both of 
the homoclinic orbits from the system. All of the periodic orbits in a neighborhood of the 
double homoclinic orbit have a stability that is determined by whether the Jacobian of the 
non-resonant saddle at the origin has a positive trace (unstable) or negative trace (stable). 
The boundaries between the parameter regions with different types of periodic orbits will 
be given by parameter values at which homoclinic orbits occur. There are four types of 
homoclinic orbits. In terms of the return map these correspond to fixed points at of ^n, 
022, (^12 ° 4>2i and 021 o 012. The parameter curves yielding the first two types of homoclinic 
orbits are simply 5i = and &2 = 0. The parameter values yielding the more complicated 
homoclinic orbits are given by 0261 + 62 = and 0162 + 61 = 0. Since the return maps are 
orientation preserving, a > 0, and taking the domain and ranges into account, these curves 
of bifurcating homoclinic orbits occur in the quadrants of the (61, 62) plane in which bi and 
62 have opposite signs. This completes the description of the unfolding of this codimension 
two bifurcation. 

In higher dimensions, double homoclinic orbits can be more complicated than the ones 
that occur for planar vector fields. The additional complication is due to the fact that 
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the return maps along two separatrices of a one-dimensional unstable manifold need not 
match in the same fashion that they do for a planar vector field. If the unstable eigenvalue 
is smaller in magnitude than all of the stable eigenvalues, then the periodic orbits that 
bifurcate from the homoclinic cycles are all stable, and there are only a finite number of 
them for any given parameter value. Indeed there are at most two. The situation that 
we describe was analyzed by Gambaudo, Glendinning and Tresser |45], who called this the 
gluing bifurcation. Their analysis is based upon the study of return maps for the double 
cycle. Let us describe a bit more. 

Let X be a vector field on i?" with an equilibrium point at the origin having a single 
unstable eigenvalue and a stable eigenspace of dimension n—1 whose spectrum lies to the 
left of the line Re{X) < —Xu- Assume further that both unstable separatrices of the origin 
are homoclinic, and that X is embedded in a generic two parameter family. Form a cross 
section S to the stable manifold Ty^(O) near the origin and let the components of 'S — W^^O) 
be Si and E2. The return maps 4>i and (j)2 of Si and S2 will be continuous, contracting and 
have images that are punctured neighborhoods of the (first) intersection of each unstable 
separatrix with S. Thus we can abstract the situation that we encounter to the study of the 
"quasicontractions" (/)i : Si ^ S and 02 : S2 — > S. If the images of these maps intersect 
VF^(O), then there will be more components formed from subsequent compositions of 4>i 
and 4'2- Using "symbolic dynamics", one can analyze the periodic orbits that can form 
from the iterated compositions of a pair of maps. The following one-dimensional model is 
adequate to represent the general situation. Consider a map / : [—1,1] [~li 1] that is 
discontinuous at the origin and contracting. The "kneading theory" or symbolic dynamics 
of / characterize its behavior. This gives a combinatorial procedure for determining the 
dynamics of / from the trajectories of the points — 1, 0~, O"*", 1. For maps of the form of 
/, there are zero, one or two periodic orbits. The patterns of signs of along periodic 

orbits are greatly restricted and can be assigned a signature that is a rational number. 
When there are two periodic orbits, their signatures p/q,p'/q' are Farey neighbors; i.e., 
\pq' — qp'\ = 1. The different families of periodic orbits appear and disappear in complex 
homoclinic orbits. Thus, even without the occurrence of horseshoes, the gluing bifurcation 
has an unfolding with an infinite number of curves in its bifurcation set. 

Other connecting orbit bifurcations have been studied by Chow, Deng and Fiedler |10|, 
Deng |13], Glendinning and Sparrow |30|, Glendinning and Tresser |31|, and Schecter ^ 



4.8 Computations of Connecting Orbits 

The numerical computation of connecting orbits has only recently been investigated, and the 
construction of robust algorithms that effectively compute connecting orbits in large classes 
of vector fields remains a challenge. We make a few comments concerning the numerical 
difficulties and fewer suggestions for how these difficulties might be confronted. We seek 
to solve the following problems. In generic one parameter families of vector fields, there 
are isolated points with connecting orbits between equilibrium points. We want accurate 
calculations of these parameter values and of the resulting orbits. (There may also be 
accumulation points of connecting orbits of increasing length and complicated topological 
structure in generic one parameter families @.) In generic two parameter families of vector 
fields, there are curves of parameter values at which connecting orbits occur. We seek to use 
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continuation methods to compute these curves. The endpoints of these curves of homochnic 
and heteroclinic bifurcations are frequently global codimension two bifurcations. We seek to 
classify these codimension two bifurcations and construct algorithms for their computation. 
The unfoldings of some of these codimension two bifurcations have been described above, but 
there is a long list of cases that have yet to be analyzed fully - even at the level of Silnikov 
approximations. Thus the goal of implementing the computation of unfoldings of global 
codimension two bifurcations involves mathematical as well as computational questions. 

Generic connecting orbits are intersections of stable and unstable manifolds of equilibria. 
Therefore, algorithms for the reliable calculation of these invariant manifolds might seem 
to form the basis for computation of the connecting orbits. This strategy is clearly an 
expensive one from a computational point of view, so we seek more direct methods for 
computing connecting orbits. In practice, most examples in higher dimensional vector fields 
have been computed by tracking periodic orbits that bifurcate from homoclinic orbits, with 
high period periodic orbits used as approximations to the homoclinic orbits. When the 
periodic orbits involved are not stable, then tracking the periodic orbits requires the use 
of algorithms for solving boundary value problems or algorithms that find fixed points of 
a cross section. Since the periodic orbits usually disappear at the homoclinic orbit, these 
computations are hard to implement in an automatic fashion. We seek methods that are 
more direct. 

Mathematically, a connecting orbit is the solution of a boundary value problem on an 
infinite interval of time. To construct algorithms based upon boundary value solvers for 
finding connecting orbits, one wants to convert the problem into one involving a finite time 
interval. This can be done approximately by introducing linear (or polynomial) approxi- 
mations for the local stable and unstable manifolds containing the ends of the connecting 
orbits and seeking trajectories that begin on the local unstable manifold and end on the local 
stable manifold. Beyn Q, Chow and Lin [11|, Doedel and Friedman |21, 22 1 and Schecter 
| ]59| have all considered algorithmic aspects of the computation of homoclinic orbits, but 
the only computations that have been examined carefully involve planar vector fields. 

One of the pragmatic questions concerning the computation of homoclinic bifurcations is 
whether the goal is to obtain an accurate approximation of the parameter values at which 
the bifurcation occurs or whether one is primarily interested in the computation of an 
accurate approximation to the homoclinic orbit. These are substantially different questions 
for reasons that we now describe. Suppose a small error has been made in the determination 
of the parameter value fj, for which there is a homoclinic orbit in a system. We ask whether 
there is a trajectory for parameter value that closely approximates the homoclinic orbit 
that exists for parameter value fj,. The distance of closest approach between the stable and 
unstable manifolds will be of order n — fi' since compact portions of the manifolds vary 
smoothly with the parameters. Thus we want to know how close the trajectory through 
a point close to the stable manifold comes to the equilibrium. Estimates can be obtained 
from linear vector fields. 

Consider a linear vector field X with a hyperbolic equilibrium point at the origin and 
assume that the stable and unstable manifolds are coordinate subspaces. The vector field 
"separates" into a stable system and an unstable system. Along trajectories, the unstable 
coordinates grow at an exponential rate. Now look at a generic point x on the stable man- 
ifold and estimate the distance of closest approach to the equilibrium point of a trajectory 
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passing through a point x' close to x. If the difference x' — x has a non-zero component in 
the direction of the strongest unstable eigenvalue of X, then this component will grow ex- 
ponentially at a rate which is the eigenvalue of the strongest unstable eigenvalue. As we saw 
above, a generic point on the stable manifold approaches the equilibrium in the direction of 
the principal (weakest) stable eigenvalue. Putting these observations together, we estimate 
the point of closest approach to the origin of the trajectory through x' as having order 
where c is the component of x' in the strongest unstable direction and /3 is the ratio of the 
magnitude of the principal stable eigenvalue to the strongest unstable eigenvalue. When 
1/P is large, then small errors in locating a point on the unstable eigenvalue lead to much 
larger distances between the closest trajectory and the homoclinic orbit that is sought. If 
we formulate a boundary value problem to find the homoclinic orbit, then we will need to 
start with a very good approximation to the homoclinic orbit to have hope of being able to 
use the boundary value solver. Trajectories depend very sensitively upon the parameters 
in a neighborhood of the equilibrium. These difficulties become still more extreme if one 
seeks the computation of a unifocal homoclinic orbit. 

We have attempted to point out some of the difficulties associated with the computation 
of connecting orbits. The computation of these dynamical structures is important due to 
their role in organizing the dynamics of a system. Consequently, this is an interesting 
problem. 



5 An Example: The Hodgkin and Huxley Equations 

This section describes the analysis of a moderately complicated dynamical system and is 
work done in collaboration with Isabel Labouriau. The Hodgkin and Huxley (HH) equations 
are a "simple" neuron model developed from experiments performed with a squid These 
equations relate the difference of electric potential across the cell membrane (V) and gating 
variables (m, n, and h) for ion channels, to the stimulus intensity (/), and temperature (T), 
as follows: 

r V = -G{V,m,n,h)+I 

I m = <^{T)[{l-m)a^{V)-mp^{V)] 

n = ^T)[{l-n)an{V)-nl3n{V)] ^'^'^> 
^ h = ^{T)[{l-h)ah{V)-hh{V)] 

where x stands for dx/dt and ^ is given by $(T) = 3(^-6.3)710 rpj^g other functions involved 
are: 

G{V, m, n, h) = m.m^h{V - Vka) + gKn\V - Vk) + g^iV - Fl) 
and the equations modeling the variation of membrane permeability: 



OLm{y) = 




draiV) = 




an{V) = 


0.1*(^) 


PniV) = 


0.125e^/80 


ah{V) = 


0.07e^/20 


l3h{V) = 


(^l + e(v+30)/io^ 


with 


^(x) = ■ 


' x/{e- - 1) 


if X / 




if x = ■ 
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Notice that ay{V) + /SyiV) ^ for all V and for y — 77X, n or h. Tlie parameters cuonj 
representing maximum conductance and equilibrium potential for the ion were obtained 
from experimental data by Hodgkin and Huxley, with the values given below: 

5Na = 120 mS/cm^ = 36 mS/cm^ gi = 0.3 mS/cm^ 
Vt^^ = -115mV Vk = 12mV VL = 10.599 mV 



The values ^Na and Vk can be controlled experimentally [42|. For the results in this pa- 



per, we use the temperature T = 6.3°C and, except where stated explicitly, all the other 
parameters involved in the HH equations have the values quoted above that we call the HH 
values. 

We describe some of the bifurcations of the HH equations as an illustration of the theory 
and algorithms described above. The local bifurcations of equilibria were calculated in terms 
of the the derivatives of the HH equations at the equilibrium points. Global bifurcations 
could only be studied by numerically integrating the HH equations. Our analysis of local 
bifurcations used the symbolic computer program Maple to implement the calculation of 
saddle-node and Hopf bifurcation curves. For y = m,n, or h the equation for y in (HH) 
is linear in y, so the last three components of an equilibrium solution (V^,, M^,, N^,, H^,) of 
(HH) can be written as functions of K : 

y*=yoo K) = — ^ , ^ for y = m,n,h . 

Substituting for y = m,n,h in the first equation, we get: 

G(K,moo(K),noo(V;),/ioo(K)) = /(K) = / . 

Thus, for fixed Vk there is exactly one value of / for which (K, m,,,, n*, h^) is at equilibrium. 
Note that derivatives of (HH) are independent of I. 

When Vk has the HH value of 12 mV, / is monotonic and (HH) has a unique equilibrium 
for each value of /. For fixed lower values of Vk, there are two saddle- node bifurcations 
as / is varied, creating a region with three equilibria. The two curves of saddle-nodes 



terminate at a cusp point. See also Q. The saddle-node curves in the I x Vk plane were 
computed parametrically with V^, as the independent parameter. The equations describing 
the saddle-node curves involve the determinant of the matrix of first derivatives of (HH) at 
an equilibrium point. We calculated an explicit expression for this determinant symbolically 
with Maple. By solving the equation that the determinant vanishes for Vk at equilibrium 
values of (V*, m*, n^,, /i,.), we obtained the curve Vk(K) of parameter values corresponding 
to zero eigenvalues. 

To determine the parameter values at which Hopf bifurcation occurs, it is necessary 
to compute eigenvalues of the matrix of first derivatives of (HH) at an equilibrium point. 
There is a pair of purely imaginary eigenvalue when the characteristic polynomial + 
csx^ + C2X^ ^ c\x ^ cq of this matrix satisfies simultaneously the third degree equation 
cf — C1C2C3 -|- C0C3 = and the inequality C1C3 > 0. These are the expressions that result 
from the Sylvester resultant calculation described earlier. Again, we computed this equation 
symbolically, assuming a given value of Kk, and solved for Vk- The graph we obtained for 
the solution of this equation and inequality disagrees slightly with the findings of Holden, et 
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al, 1 44 1 for the HH value = 36. Takens-Bogdanov bifurcations occur when the equations 
defining Hopf bifurcations and saddle-node bifurcations are satisfied simultaneously. 

The saddle-node and Hopf bifurcations are the only codimension one bifurcations that 
can be computed explicitly from (HH) without numerical integration. The presence of dou- 
ble cycles where the two periodic orbits created at Hopf bifurcation points coalesce and 
disappear has been established previously |3^, ^ and the existence of saddle loops 
emanating from the Takens-Bogdanov points is predicted by bifurcation theory [^6|. To 
determine further information about global bifurcations, we rely upon numerical integra- 
tions that were performed with the computer program DsTool This program establishes 
a graphical interface and display for investigating bifurcations of dynamical systems. It 
allows one to mark points in a two-dimensional parameter space with identifying symbols 
and to display phase portraits that correspond to these points. The computed data for local 
bifurcations was displayed in the parameter space window. By searching for the boundaries 
of parameter regions yielding structurally stable dynamics and using our knowledge of the 
unfoldings of codimension two bifurcations, we deduced the location of curves along which 
global bifurcations take place. We obtain a consistent picture of the bifurcation diagrams 
for (HH) in the two-dimensional / x Vk parameter plane. These diagrams have not been 
proved to be correct, but they are based upon strong numerical evidence. 

The bifurcation diagram resulting from our numerical investigations is shown in Figure |^. 
Its main features include a curve of double cycles (dc) which enters the cusp region with 
three equilibrium points and terminates at a degenerate Hopf bifurcation (dh) close to 
the Takens-Bogdanov point (tb). These double cycles are the ones described in |Q. The 
curve of saddle loops (si) emanating from the Takens-Bogdanov point crosses the Hopf 
curve beyond the degenerate Hopf point, and then turns sharply. From this sharp bend, 
it proceeds almost parallel to the saddle-node curve (sn). The saddle loops appear to 
undergo a reversal of orientation along this portion of the curve (si). After the orientation 
has reversed, one encounters a set of parameters at which the unstable eigenvalue and the 
weakest stable eigenvalue have equal magnitudes. This point (tsl) in parameter space is a 
twisted neutral saddle loop, and there are additional curves of untwisted twice rounding 
and period doubling bifurcations that emerge from (tsl). 

When is decreased from the HH value of 36 mS/cm^ the Takens-Bogdanov point 
in the I xVk plane moves towards the cusp point and past it. This agrees qualitatively 
with the findings of [0], but their results differ from ours in the value of for which 
the Takens-Bogdanov point moves past the cusp. The unfolding for the codimension three 
bifurcation in which cusp and Takens-Bogdanov bifurcations coincide has been analyzed by 
Guckenheimer |34] and Dumortier and Roussarie |16|. The geometry of the unfolding of 
this codimension three bifurcation can be visualized by drawing a two dimensional sphere 
that encloses the codimension three point in the three dimensional parameter space of the 
unfolding |34]. 

To further explore the effect of this codimension three bifurcation on the bifurcation 
diagrams of (HH), we decreased from the HH value of 36 mS/cm^ to 12 mS/cm^ and 
computed another bifurcation diagram in the / x Vk plane. The new bifurcation diagram is 
shown in Figure Among its features are a curve of double cycles (dc) that terminates at 
a neutral saddle loop point (nsl) instead of a double Hopf bifurcation as in the unfolding of 
the codimension three bifurcation. The point (nsl) does not lie on the saddle loop branch 
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emanating from the Takens-Bogdanov point (tb), however. Instead it ends on a saddle loop 
that encloses both equilibrium points. This branch of saddle loops ends on both branches 
of saddle- nodes at saddle node loops (snl). The branch of twisted saddle loops (tsl) that 
was present at the higher value of remains. It starts on the saddle-node curve at another 
saddle node loop. The twisted saddle loops still passes through a neutral point (tnsl) at 
which bifurcation curves of period doublings (pd) and doubled saddle loops (dsl) originate. 
Our proposed bifurcation diagrams for (HH) near the cusp points appear to be compatible 
with the unfolding of the Takens-Bogdanov cusp codimension three bifurcation, though 
the diagrams drawn here are sufficiently far from the codimension three bifurcation that 
significant differences with its unfolding exist. 

The intent of this rapid tour of the complicated dynamics of a realistic biological model 
was to draw attention to the interplay between theory and numerical exploration. Since 
many regions in the parameter space are very small, the theory helps to guide the exploration 
to discover the correct bifurcation diagram. 

6 The Effects of Symmetry 

In the previous sections, non-generic bifurcation has arisen in multiparameter families of 
vector fields and is a codimension two, or higher, phenomenon. Non-generic bifurcation is 
often encountered in the real world as a result of symmetry and in this section we discuss 
the bifurcation of symmetric vector fields. 

Within a generic one parameter family of vector fields, x = fx{x), bifurcations of equi- 
libria are either saddle-node bifurcations, when a single eigenvalue passes through 0, or 
Hopf bifurcations, when a pair of complex conjugate eigenvalues cross the imaginary axis. 
These are both well understood [^6|, but with symmetry, generic bifurcation can be much 
more complicated. One of our goals has been to explore how complex local bifurcation with 
symmetry can be. 

We first make a few definitions. Let G be a finite subgroup of 0(n) acting on i?". An 
equation x = f{x) is G-equivariant or symmetric if 

figx) = go fix) for all g £ G . 

For a fixed group G, we wish to study generic bifurcation within the class of G-equivariant 
vector fields. 

We now make two assumptions on the group action. The first is that the action of G is 
absolutely irreducible. This means that the only matrices which commute with G are scalar 
multiples of the identity. The second is that the action of G is fixed point free: the only 
point fixed by G is the origin. These assumptions have two consequences for G-equivariant 
vector fields. The first is that the origin is always an equilibrium point of a G-equivariant 
vector field, since 

/(O) = figO) = 5/(0) => /(O) = . 

The second consequence is that the linearization at the origin is a multiple of the identity. 
This follows from the calculations: 



D{fog)\o=Df\g^o)9 = Df\og 
and Dif o g)\o = Dig f)\o=gDf\o 
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Because -D/|o commutes with all g £ G and G is absolutely irreducible, D/|o must be a 
scalar multiple of the identity. We may now state the following: 

Theorem 6.1 Let G be fixed point free and absolutely irreducible and let x = f\{x) be a one- 
parameter family of G-equivariant vector fields. If the origin undergoes a bifurcation with 
one zero eigenvalue for A = Aq, then the linearization at the origin is given by -D/a|(o,Ao) ~ 0- 
In other words, all the eigenvalues must pass through zero simultaneously. 

This degenerate linearization occurs as a codimension one phenomenon due to the symmetry 
of the problem. We wish to analyze the bifurcations at the origin, but the possibilities are 
not clear and there are still many restrictions due to the equivariance. 

To study the bifurcation in an analogous fashion as for the non-symmetric case, we will 
compute a normal form by computing a Taylor series at the origin: 

oo 

/a(x) = ^P,(x). 

i=l 

Each Pi is a homogeneous, degree i, G-equivariant. The Pi are not always easy to calculate 
and we are led to the study of polynomial invariants and equivariants which combines 



elements of combinatorial theory and modern commutative algebra |^, 64], but we shall 
not delve into this topic in this paper. 

We shall now study a simple example: -D4 acting on R^, which may be thought of as 
the symmetry of the square in the plane. Generators for this action are given by 

and 





and it is easy to show that this action is fixed point free and absolutely irreducible. We 
look at the degree three power series expansion at the origin: 



X 



Xx + ax^ + bxy^ 
y = Xy + ay^ + bx^y 

As A passes through 0, the origin changes stability, but what else happens in the process? 
We can get some additional information from looking at subgroups of G and fixed point 
subspaces. 

li H <Z G, then we define the fixed point subspace of the subgroup H by 

Fix(i7) = {x e i?" such that hx = x for ah he H}. 

Fixed point subspace are important because they are invariant under the flow. This is easy 
to see because if x S Fix{II), then f{x) = f{hx) = hf{x) so f{x) G Fix(if), i.e., the vector 
field on Fix(ff) is tangent to Fix(i/). The fixed point subspaces for the example are the 
axes of symmetry plus the trivial ones: the origin and all of E?. 

To explain bifurcation in some of the fixed point subspaces, Cicogna and Vander- 
bauwhede formulated the following lemma: 
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Lemma 6.2 (The Fixed Point Branching Lemma) jp^l If V is a one- dimensional fixed 
point subspace, then there is a branch of solutions in V in i?" x R that pass through (0, 0) 
and is part of a generic bifurcation. 

For the D4 example, this means that on the axes of symmetry there are pitchfork bifur- 
cations. For this simple problem, we may verify this result by explicitly computing the 
equilibrium points which are given by: 

(0,0) 



-A 



-A 



a + \ a + b 



So if a < and a + 6 < 0, then as A passes through zero, eight new equilibrium points are 
created. It is typical behavior to get multiple equilibrium points in equivariant bifurcations 
and a substantial theory has been developed by M. Field |17|, M. Golubitsky |33] and 
others. This work, based on group theory and singularity theory, can give information 
about bifurcation to equilibrium points, period orbits and quasiperiodic orbits. However, 
this is a far from complete categorization of the behavior which may be obtained from more 
complicated group actions. 

Now consider the action of a group G on i?^ generated by 



/ 


-1 









/ 





1 









1 


;] 


and 












V 












I 


1 








The action is fixed point free, absolutely irreducible and, a few short calculations show, the 
degree three equivariant vector field has the form: 



Xl 

X3 



(A + aixf + 02X2 + a3x|)xi 
(A + 01X2 + a2xl + a3xl)x2 
(A + 01X3 + a2xl + a3xl)x3 



(2) 



One of the consequences of the symmetry is the invariance of both the coordinate axes and 
the coordinate planes. Guckenheimer and Holmes [35| observed that this system of equations 
has structurally stable heteroclinic cycles for open regions of the parameter space. 

Theorem 6.3 // either 03 < ai < 02 < or 02 < ai < 03 < and 2ai > 02 + 03 then 
Equations |^ have attracting heteroclinic cycles for A > and a globally attracting equilibrium 
at the origin for A < 0. 

This theorem essentially follows from examining the flow in the invariant coordinate planes. 
The consequence is that within families of G-equivariant vector fields there are generic one 
parameter bifurcations from an attracting equilibrium point to structurally stable attracting 
limit cycles. We continue to study more complicated symmetry groups in an attempt to 
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discover how complicated post bifurcation behavior can be. The answer is that a stable 



equilibrium point can bifurcate directly to a chaotic attractor of small amplitude |38] 
Consider the action of G on R'^ generated by 
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-1 








M 







1 
















1 





V 














and 



/ 1 \ 

10 

1 

\ 1 / 



as studied by Field and Swift |19]. The cubic truncation of the general G-equivariant vector 
field is given by 

Xl = (A + aixf + a2X2 + 03X3 + a4x|)xi 
X2 = (A + aixl + 02X3 + a^xl + a4x'l)x2 
X3 = (A + aixl + a2xl + a^xj + 04X2)^3 
X4 = (A + 01X4 + a2X^ + 03X2 + a4X3)x4 

First observe that if all the Oj = —1, then in the flow of the vector field, for A < the origin 
is globally attracting, and for A > there is an invariant attracting sphere of radius \/A. 
The sphere is normally hyperbolic. If the Oj are close to —1, then there will be an invariant 
topological 3-sphere in the post bifurcation flow (The Invariant Sphere Theorem [|l^]). 

We may now ask about the dynamics on the invariant 3-sphere. We rescale the phase 
space variables by \/A and the independent variable by 1/A to blow up the sphere, which 
is equivalent to setting A = 1 in the above equations. Varying the parameters Oj, we study 
the flow on the sphere realizing that all the dynamics on the sphere represent possible post 
bifurcation behavior of the vector field. Structurally stable objects give persistent post 
bifurcation behavior for G-equivariant families. All indications imply that this system has 
no chaotic behavior, but consider G C G, the elements of the group G which are orientation 
preserving. A new term must be added to the equivariant vector field which results in 



Xl 
X2 

X3 
X4 



(A -|- aixf + 02X2 + a3x| -|- a4X4)xi -|- 6X2X3X4 
(A -|- 01X2 + a2X^ + 03X4 -|- a4xf)x2 — 6X1X3X4 
(A -|- aix| -|- 02X4 -|- a^xf + a4X2)x3 -|- 6X1X2X4 
(A -|- 01X4 -|- a2xf -|- 03X2 -|- a4x|)x4 — 6X1X2X3 



The geometry associated with this group action appears to be very rich and was first 
examined by Guckenheimer and Worfolk [^8|. For Oj ~ — 1 and A > 0, there is an invariant, 
attracting, topological sphere in the flow on which we wish to study the dynamics. Also, 
as in the earlier examples, there are structurally stable heteroclinic cycles given by the 
intersection of the invariant sphere and the invariant Xj — Xj+i orthogonal coordinate planes. 

The geometry associated with this group action appears to be very rich. We may think of 
objects as bifurcating from the heteroclinic cycles as they are broken by the introduction of 
fixed points in the cycles. A numerical study reveals unmistakable signs of chaotic behavior: 
period doubling cascades and Silnikov homoclinic orbits. This leads to the conclusion that 
the flow on the invariant sphere can contain chaotic attractors. All hyperbolic invariant 
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sets are possible post bifurcation behavior, hence we can expect to bifurcate directly from a 
trivial equilibrium to a chaotic attractor of small amplitude. We are currently applying the 
techniques presented for the analysis of homoclinic bifurcations to the task of analytically 
proving the existence of chaotic behavior in a neighborhood of the heteroclinic cycles. This 
would verify the following conjecture formulated from the numerical exploration. 

Conjecture 6.4 There is an open regionU in the space of one parameter families of vector 
fields on equivariant with respect to G, such that X\ G U implies that a bifurcation of 
the trivial equilibrium point occurs in Xx at \ = \q that produces "instant chaos". This 
means that if U is a neighborhood of the origin in i?^, then there is an e > such that Xx 
has a chaotic hyperbolic invariant set contained in U for < X — Xq < e. 
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8 Captions for Figures 

Figure |l|: Evolve the stable manifold of p using geodesic curves. 

Figure |2|: The stable manifold of the origin in the Lorenz system projected onto the x-z 
plane. 

Figure |3|: Depiction of section of W^~^ H W^~^ for typical homoclinic orbits. 
Figure |^: Unifocal and bifocal homoclinic orbits. 

Figure |5|: Bifurcation diagram for an untwisted resonant homoclinic orbit. 
Figure |6|: Bifurcation diagrams for twisted resonant homoclinic orbit. 
Figure |^: Return maps for the planar gluing bifurcation. 
Figure |8|: Bifurcation diagram for the Hodgkin and Huxley equations. 



Figure^ Bifurcation diagram for the Hodgkin and Huxley equations, with gx = 12 mS/cm^. 



30 



J. Guckenheimer and P. Worfolk 



References 

[1] F. Argoul, A. Arneodo, and P. Richetti. Experimental evidence for homoclinic chaos 
in the Belousov-Zhabotinskii reaction. Phys. Lett. A, 120(6) :269-275, 1987. 

[2] V. I. Arnold. Geometrical Methods in the Theory of Ordinary Differential Equations. 
Springer- Verlag, New York, 1988. 

[3] A. Back, J. Guckenheimer, M. Myers, F. Wicklin, and P. Worfolk. dstool: Computer 
assisted exploration of dynamical systems. Notices Amer. Math. Soc, 39(4):303-309, 
1992. 

[4] C. Baesens, J. Guckenheimer, S. Kim, and R. Mackay. Three coupled oscillators: 
Mode-locking, global bifurcations and toroidal chaos. Phys. D, 49:387-475, 1991. 

[5] L. Belyakov. Bifurcation of systems with homoclinic curve of a saddle-focus with saddle 
quantity zero. Mat. Zametki, 36(5):681-689, 1985. 

[6] M. Benedicks and L. Carleson. The dynamics of the Henon map. Ann. of Math., 
133:73-169, 1991. 

[7] W. J. Beyn. The numerical computation of connecting orbits in dynamical systems. 
IMA J. Numer. Anal, 9:169-181, 1990. 

[8] M. Bosch and C. Simo. Attractors in a Silnikov-Hopf scenario and a related one- 
dimensional map. Phys. D, 62:217-229, 1993. 

[9] S.-N. Chow, B. Deng, and B. Fiedler. Homoclinic bifurcation at resonant eigenvalues. 
J. Dynamics Differential Equations, 2(2):177-244, 1990. 

[10] S.-N. Chow, B. Deng, and D. Terman. The bifurcation of homoclinic and periodic 
orbits from two heteroclinic cycles. SIAM J. Math. Anal, 21(l):179-204, 1990. 

[11] S. N. Chow and X. B. Lin. Bifurcation of a homoclinic orbit with a saddle-node 
equilibrium. Differential Integral Equations, 3:435-466, 1990. 

[12] B. Deng. Homoclinic bifurcations with nonhyperbolic equilibria. SIAM J. Math. Anal, 
21(3):693-720, 1990. 

[13] B. Deng. The bifurcations of countable connections from a twisted homoclinic loop. 
SIAM J. Math. Anal, 22(3):653-679, 1991. 

[14] B. Deng. Homoclinic twisting bifurcation and cusp horseshoe maps. Preprint, Dec. 
1991. 

[15] M. H. Dulac. Sur les cycles limites. Bull Soc. Math. Anal, 51:45-188, 1923. 

[16] F. Dumorticr, R. Roussaric, and J. Sotomayor. Generic 3-paramctcr families of vector 
fields on the plane, unfolding a singularity with nilpotcnt linear part. The cusp case of 
codimension 3. Ergodic Theory Dynamical Systems, 7:375-413, 1987. 



Dynamical Systems: Some Computational Problems 



31 



[17] M. Field. Equivariant dynamics. Contemp. Math., 56:69-96, 1986. 

[18] M. Field. Equivariant bifurcation theory and symmetry breaking. J. Dynamics Dif- 
ferential Equations, 1(4):369-421, 1989. 

[19] M. Field and J. Swift. Stationary bifurcation to limit cycles and heteroclinic cycles. 
Nonlinearity, 4:1001-1043, 1991. 

[20] A. C. Fowler and C. T. Sparrow. Bifocal homoclinic orbits in four dimensions. Non- 
linearity, 4:1159-1182, 1991. 

[21] M. J. Friedman. Numerical analysis and accurate computation of heteroclinic orbits in 
the case of center manifolds. To appear in J. Dynamics Differential Equations, 1992. 

[22] M. J. Friedman and E. J. Doedel. Computational methods for global analysis of ho- 
moclinic and heteroclinic orbits: a case study. To appear in J. Dynamics Differential 
Equations, 1992. 

[23] A. T. Fuller. Conditions for a matrix to have only characteristic roots with negative 
real parts. J. Math. Anal. AppL, 23:71-98, 1968. 

[24] P. Gaspard. Generation of a countable set of homoclinic flows through bifurcation. 
Phys. Lett. A, 97(l):l-4, 1983. 

[25] P. Gaspard, R. Kapral, and G. Nicolis. Bifurcation phenomena near homoclinic sys- 
tems: A two-parameter analysis. J. Statist. Phys., 35(5/6) :697-727, 1984. 

[26] P. Gaspard and X.-J. Wang. Homoclinic orbits and mixed-mode oscillations in far- 
from-equilibrium systems. J. Statist. Phys., 48(1):151-199, 1987. 

[27] P. Glendinning. Travelling wave solutions near isolated double-pulse solitary waves of 
nerve axon equations. Phys. Lett. A, 121:411-413, 1987. 

[28] P. Glendinning. Subsidiary bifurcations near bifocal homoclinic orbits. Math. Proc. 
Cambridge Philos. Soc, 105:597-605, 1989. 

[29] P. Glendinning and C. Sparrow. Local and global behavior near homoclinic orbits. J. 
Stat. Phys., 34(5/6):645-696, 1984. 

[30] P. Glendinning and C. Sparrow. T-points: A codimension two heteroclinic bifurcation. 
J. Stat. Phys., 43(3/4) :479-488, 1986. 

[31] P. Glendinning and C. Tresser. Heteroclinic loops leading to hyperchaos. J. Physique 
Lett., 46(8):347-352, 1985. 

[32] M. Golubitsky and V. Guillemin. Stable Mappings and Their Singularities. Springer- 
Verlag, New York, 1973. 

[33] M. Golubitsky, I. Stewart, and D. Shaeffer. Singularities and Groups in Bifurcation 
Theory, volume II. Springer- Verlag, New York, 1988. 



32 J. Guckenheimer and P. Worfolk 

[34] J. Guckenheimer. Multiple bifurcation problems for chemical reactors. Phys. D, 20:1- 
20, 1986. 

[35] J. Guckenheimer and P. Holmes. Structurally stable heteroclinic cycles. Math. Proc. 
Cambridge Philos. Soc, 103:189-192, 1988. 

[36] J. Guckenheimer and P. J. Holmes. Nonlinear Oscillations, Dynamical Systems, and 
Bifurcations of Vector Fields. Springer- Verlag, New York, 1983. 

[37] J. Guckenheimer and R. Williams. Structural stability of Lorenz attractors. Publ. 
Math. IHES, 50:59-72, 1979. 

[38] J. Guckenheimer and P. Worfolk. Instant chaos. Nonlinearity, 5:1211-1222, 1992. 

[39] B. D. Hassard and L.-J. Shiau. Isolated periodic solutions of the Hodgkin-Huxley 
equations. J. Theoret. Biol, 136:267-280, 1989. 

[40] M. Henon. A two-dimensional mapping with a strange attractor. Comm. Math. Phys., 
50:69-77, 1976. 

[41] P. Hirschberg and E. Knobloch. Silnikov-Hopf bifurcation. Phys. D, 62:202-216, 1993. 

[42] A. L. Hodgkin and A. F. Huxley. Current carried by sodium and potassium ions 
through the membrane of the giant axon of loligo. J. Physiology, 116:449-472, 1952. 

[43] A. L. Hodgkin and A. F. Huxley. A quantitative description of membrane current 
and its applications to conduction and excitation in nerve. J. Physiology, 117:500-544, 
1952. 

[44] A. V. Holden, M. A. Muhamad, and A. K. Schierwagcn. Repolarizing currents and 
periodic activity in nerve membrane. J. Theoret. NeurobioL, 4:61-71, 1985. 

[45] P. Glendinning J. M. Gambaudo and C. Tresser. Gollage de cycles et suites de Farey. 
C. R. Acad. Sci. Paris, 299:711-714, 1984. 

[46] H. B. Keller. Numerical solution of bifurcation and nonlinear eigenvalue problems. In 
P. Rabinowitz, editor. Applications of Bifurcations Theory, pages 359-384. Academic 
Press, 1977. 

[47] H. Kokubu. Homoclinic and heteroclinic bifurcations of vector fields. Japan J. Appl. 
Math, 5:455-501, 1988. 

[48] I. S. Labouriau. Degenerate Hopf bifurcation and nerve impulse. Part II. SIAM J. 
Math. Anal, 20:1-12, 1989. 

[49] R. Loos. Generalized polynomial remainder sequences. In G.E Collins B. Buchberger 
and R. Loos, editors. Computer Algebra - Symbolic and Algebraic Computation, pages 
115-137. Springer- Verlag, 1982. 

[50] E. N. Lorenz. Deterministic non-periodic flow. J. Atmospheric Sci., 20:130-141, 1963. 



Dynamical Systems: Some Computational Problems 



33 



[51] V. Luk'yanov. Bifurcations of dynamical systems with a saddle-point-separatrix loop. 
J. Differential Equations, 18:1049-1059, 1983. 

[52] L. Mora and M. Viana. Abundance of strange attractors. Preprint, 1991. 

[53] J. Palis and W. de Melo. Geometric Theory of Dynamical Systems. Springer- Ver lag. 
New York, 1982. 

[54] R. Richetti, F. Argoul, and A. Arneodo. Type-II intermittency in a periodically driven 
nonlinear oscillator. Phys. Rev. A, 34(l):726-729, 1986. 

[55] J. Rinzel and R. N. Miller. Numerical solutions of the Hodgkin-Huxley equations. 
Math. Biosc, 49:27-59, 1980. 

[56] S. Schecter. The saddle-node separatrix-loop bifurcation. SI AM J. Math. Anal., 
18(4):1142-1156, 1987. 

[57] S. Schecter. Simultaneous equilibrium and heteroclinic bifurcation of planar vector 
fields via the Melnikov integral. Nonlinearity, 3:79-99, 1990. 

[58] S. Schecter. singularity theory and heteroclinic bifurcation with a distinguished 
parameter. J. Differential Equations, 99:306-341, 1992. 

[59] S. Schecter. Numerical computation of saddle-node homoclinic bifurcation points, to 
appear in SIAM J. Numer. Anal., 1993. 

[60] L. P. Silnikov. A case of the existence of a countable number of periodic motions. 
Soviet Math. Dokl, 6:163-166, 1965. 

[61] L. P. Silnikov. The existence of a denumerable set of periodic motions in four- 
dimensional space in an extended neighborhood of a saddle-focus. Soviet Math. Dokl., 
8(l):54-58, 1967. 

[62] L. P. Silnikov. On the generation of a periodic motion from trajectories doubly asymp- 
totic to an equilibrium state of saddle type. Math. USSR-Sb., 6(3):427-438, 1968. 

[63] L. P. Silnikov. A contribution to the problem of the structure of an extended neigh- 
borhood of a rough equilibrium of saddle-focus type. Math. USSR-Sb., 10(1):91-102, 
1970. 

[64] R.P. Stanley. Invariants of finite groups and their applications to combinatorics. Bull. 
Amer. Math. Soc, 1(3):475-511, 1979. 

[65] D. Terman. The transition from bursting to continuous spiking in excitable membrane 
models. J. Nonlinear Science, 2(2):135-182, 1992. 

[66] C. Tresser. About some theorems by L.P. Silnikov. Ann. Inst. H. Poincare Sect. A, 
40(4):441-461, 1984. 




gure 2: The stable manifold of the origin in the Lorenz system projected onto the x-z 
ane. 
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Figure 3: Depiction of section of W^'^ fl 1^"+ for typical homoclinic orbits 
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Figure 4: Unifocal and bifocal homoclinic orbits. 
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Figure 5: Bifurcation diagram for an untwisted resonant homoclinic orbit. 



twice rounding 
homoclinic orbit 



1 periodic orbit ^ 



7 



homoclinic orbit 




/ 2 periodic orbits 

\ \ 1 periodic orbit 

A 

period doubling bifurcation 



Figure 6: Bifurcation diagrams for twisted resonant homoclinic orbit. 
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Figure 7: Return maps for the planar gluing bifurcation. 



4 





Figure 8: Bifurcation diagram for the Hodgkin and Huxley equations. 



5 




6 



